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Abstract 

HYDJET++ is a Monte-Carlo event generator for simulation of relativistic heavy 
ion AA collisions considered as a superposition of the soft, hydro-type state and the 
hard state resulting from multi-parton fragmentation. This model is the develop- 
ment and continuation of HYDJET event generator (Lokhtin Sz Snigirev, 2006, 
EPJC, 45, 211). The main program is written in the object-oriented C++ language 
under the ROOT environment. The hard part of HYDJET++ is identical to the 
hard part of Fortran-written HYDJET and it is included in the generator struc- 
ture as a separate directory. The soft part of HYDJET++ event is the "thermal" 
hadronic state generated on the chemical and thermal freeze-out hypersurfaces ob- 
tained from the parameterization of relativistic hydrodynamics with preset freeze- 
out conditions. It includes the longitudinal, radial and elliptic flow effects and the 
decays of hadronic resonances. The corresponding fast Monte-Carlo simulation pro- 
cedure, C++ code FAST MC (Amelin et al., 2006, PRC, 74, 064901; 2008, PRC, 77, 
014903) is adapted to HYDJET++. It is designed for studying the multi-particle 
production in a wide energy range of heavy ion experimental facilities: from FAIR 
and NIC A to RHIC and LHC. 
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PROGRAM SUMMARY 



Manuscript Title: Heavy ion event generator HYDJET++ (HYDrodynamics plus 
JETs) 

Authors: Igor Lokhtin, Ludmila Malinina, Sergey Petrushanko, Alexander Snigirev, 
Ionut Arsene, Konrad Tywoniuk 

Program Title: HYDJET++, version 2.0 

Programming language: C++ (however there is a Fortran-written part which is 
included in the generator structure as a separate directory) 

Size of the package: 3.5 MBytes directory and 800 kBytes compressed distribution 
archive (without ROOT libraries). The output file created by the code in ROOT tree 
format for 100 central (0-5%) Au+Au events at y/s = 200 A GeV (Pb+Pb events 
at \fs = 5500A GeV) with default input parameters requires 40 (190) MBytes of 
the disk space. 

Computer: PC - hardware independent (both C++ and Fortran compilers and 
ROOT environment [1] should be installed) 

Operating system: Linux (Scientific Linux, Red Hat Enterprise, FEDORA, etc.) 
RAM: 50 MBytes (determined by ROOT requirements) 
Number of processors used: 1 

Keywords: Monte-Carlo event generators, relativistic heavy ion collisions, hydrody- 
namics, QCD jets, partonic energy loss, flow, quark-gluon plasma 

PACS: 24.10.Lx, 24.85.+p, 25.75.-q, 25.75.Bh, 25.75.Dw, 25.75.Ld, 25.75.Nq 

Classification: 11.2 Elementary particle physics, phase space and event simulation 

External routines /libraries: ROOT (any version) 

Subprograms used: PYTHIA event generator (version 6.401 or later), PYQUEN 
event generator (version 1.5 or later) 

Nature of the physical problem: 

The experimental and phenomenological study of multi-particle production in rela- 
tivistic heavy ion collisions is expected to provide valuable information on the dy- 
namical behaviour of strongly-interacting matter in the form of quark-gluon plasma 
(QGP) [2-4], as predicted by lattice Quantum Chromodynamics (QCD) calcula- 
tions. Ongoing and future experimental studies in a wide range of heavy ion beam 
energies require the development of new Monte-Carlo (MC) event generators and im- 
provement of existing ones. Especially for experiments at the CERN Large Hadron 
Collider (LHC), implying very high parton and hadron multiplicities, one needs 
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fast (but realistic) MC tools for heavy ion event simulations [5-7]. The main ad- 
vantage of MC technique for the simulation of high-multiplicity hadroproduction 
is that it allows a visual comparison of theory and data, including if necessary the 
detailed detector acceptances, responses and resolutions. The realistic MC event 
generator has to include maximum possible number of observable physical effects, 
which are important to determine the event topology: from the bulk properties of 
soft hadroproduction (domain of low transverse momenta pt ^ lGeV/c) such as 
collective flows, to hard multi-parton production in hot and dense QCD-matter, 
which reveals itself in the spectra of high-py particles and hadronic jets. Moreover, 
the role of hard and semi-hard particle production at LHC can be significant even 
for the bulk properties of created matter, and hard probes of QGP became clearly 
observable in various new channels [8-11]. In majority of the available MC heavy 
ion event generators, the simultaneous treatment of collective flow effects for soft 
hadroproduction and hard multi-parton in-medium production (medium-induced 
partonic rescattering and energy loss, so called "jet quenching") is lacking. Thus, in 
order to analyze existing data on low and high-p^ hadron production, test the sen- 
sitivity of physical observables at the upcoming LHC experiments (and other future 
heavy ion facilities) to the QGP formation, and study the experimental capabilities 
of constructed detectors, the development of adequate and fast MC models for si- 
multaneous collective flow and jet quenching simulations is necessary. HYDJET++ 
event generator includes detailed treatment of soft hadroproduction as well as hard 
multi-parton production, and takes into account known medium effects. 

Solution method: 

A heavy ion event in HYDJET+- 1- is a superposition of the soft, hydro-type state and 
the hard state resulting from multi-parton fragmentation. Both states are treated 
independently. HYDJET++ is the development and continuation of HYDJET MC 
model [12]. The main program is written in the object-oriented C++ language un- 
der the ROOT environment [1]. The hard part of HYDJET++ is identical to the 
hard part of Fortran- written HYDJET [13] (version 1.5) and is included in the 
generator structure as a separate directory. The routine for generation of single 
hard NN collision, generator PYQUEN [12,14], modifies the "standard" jet event 
obtained with the generator PYTHIA_6.4 [15]. The event-by-event simulation pro- 
cedure in PYQUEN includes 1 ) generation of initial parton spectra with PYTHIA 
and production vertexes at given impact parameter; 2) rescattering-by-rescattering 
simulation of the parton path in a dense zone and its radiative and collisional en- 
ergy loss; 3) final hadronization according to the Lund string model for hard partons 
and in-medium emitted gluons. Then the PYQUEN multi-jets generated according 
to the binomial distribution are included in the hard part of the event. The mean 
number of jets produced in an AA event is the product of the number of binary NN 
sub-collisions at a given impact parameter and the integral cross section of the hard 
process in NN collisions with the minimum transverse momentum transfer p™ in . In 
order to take into account the effect of nuclear shadowing on parton distribution 
functions, the impact parameter dependent parameterization obtained in the frame- 
work of Glauber- Gribov theory [16] is used. The soft part of HYDJET++ event is 
the "thermal" hadronic state generated on the chemical and thermal freeze-out hy- 
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persurfaces obtained from the parameterization of relativistic hydrodynamics with 
preset freeze-out conditions (the adapted C++ code FAST MC [17,18]). Hadron 
multiplicities are calculated using the effective thermal volume approximation and 
Poisson multiplicity distribution around its mean value, which is supposed to be 
proportional to the number of participating nucleons at a given impact parameter 
of AA collision. The fast soft hadron simulation procedure includes 1) generation 
of the 4-momentum of a hadron in the rest frame of a liquid element in accordance 
with the equilibrium distribution function; 2) generation of the spatial position of a 
liquid element and its local 4- velocity in accordance with phase space and the char- 
acter of motion of the fluid; 3) the standard von Neumann rejection/acceptance 
procedure to account for the difference between the true and generated probabil- 
ities; 4) boost of the hadron 4-momentum in the center mass frame of the event; 
5) the two- and three-body decays of resonances with branching ratios taken from 
the SHARE particle decay table [19]. The high generation speed in HYDJET++ is 
achieved due to almost 100% generation efficiency of the "soft" part because of the 
nearly uniform residual invariant weights which appear in the freeze-out momen- 
tum and coordinate simulation. Although HYDJET++ is optimized for very high 
energies of RHIC and LHC colliders (c.m.s. energies of heavy ion beams \/s = 200 
and 5500 GeV per nucleon pair respectively), in practice it can also be used for 
studying the particle production in a wider energy range down to y/s ~ 10 GeV per 
nucleon pair at other heavy ion experimental facilities. As one moves from very high 
to moderately high energies, the contribution of the hard part of the event becomes 
smaller, while the soft part turns into just a multi-parameter fit to the data. 

Restrictions: 

HYDJET++ is only applicable for symmetric AA collisions of heavy (A > 40) ions 
at high energies (c.m.s. energy \fs > 10 GeV per nucleon pair). The results obtained 
for very peripheral collisions (with the impact parameter of the order of two nucleus 
radii, b ~ 2Ra) and very forward rapidities may be not adequate. 

Running time: 

The generation of 100 central (0-5%) Au+Au events at y/s = 200A GeV (Pb+Pb 
events at y/s = 5500^4 GeV) with default input parameters takes about 7 (85) 
minutes on a PC 64 bit Intel Core Duo CPU @ 3 GHz with 8 GB of RAM memory 
under Red Hat Enterprise. 

Accessibility: 

http:/ /cern.ch/lokhtin/hydjet++ 
References for the physics model: 

[1] LP Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211. 

[2] N.S. Amelin, R. Lednicky, T.A. Pocheptsov, LP. Lokhtin, L.V. Malinina, A.M. 
Snigirev, Iu.A. Karpenko and Yu.M. Sinyukov, Phys. Rev. C 74 (2006) 064901. 

[3] N.S. Amelin, I. Arsene, L. Bravina, Iu.A. Karpenko, R. Lednicky, LP. Lokhtin, 
L.V. Malinina, A.M. Snigirev and Yu.M. Sinyukov, Phys. Rev. C 77 (2008) 
014903. 
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1 Introduction 



One of the basic tasks of modern high energy physics is the study of the fun- 
damental theory of strong interaction (Quantum Chromodynamics, QCD) in 
new, unexplored extreme regimes of super-high densities and temperatures 
through the investigation of the properties of hot and dense multi-parton and 
multi-hadron systems produced in high-energy nuclear collisions [20-22]. In- 
deed, QCD is not just a quantum field theory with an extremely rich dynami- 
cal content (such as asymptotic freedom, chiral symmetry, non-trivial vacuum 
topology, strong CP violation problem, colour superconductivity), but per- 
haps the only sector of the Standard Model, where the basic features (as phase 
diagram, phase transitions, thermalisation of fundamental fields) may be the 
subject to scrutiny in the laboratory. The experimental and phenomenological 
study of multi-particle production in ultrarelativistic heavy ion collisions [2- 
4] is expected to provide valuable information on the dynamical behaviour 
of QCD matter in the form of a quark-gluon plasma (QGP), as predicted by 
lattice calculations. 

Experimental data, obtained from the Relativistic Heavy Ion Collider (RHIC) 
at maximum beam energy in the center of mass system of colliding ions 
yHi = 200 GeV per nucleon pair, supports the picture of formation of a strongly 
interacting hot QCD matter ( "quark-gluon fluid" ) in the most central Au+ Au 
(and likely Cu+Cu) collisions [23-26]. This appears as significant modification 
of properties of multi-particle production in heavy ion collisions as compared 
with the corresponding proton-proton (or peripheral heavy ion) interactions. 
In particular, one of the important perturbative ("hard") probes of QGP is 
the medium-induced energy loss of energetic partons, so called "jet quench- 
ing" [27], which is predicted to be very different in cold nuclear matter and 
in QGP, and leads to a number of phenomena which are already seen in the 
RHIC data on the qualitative level, such as suppression of high-p-p hadron 
production, modification of azimuthal back-to-back correlations, azimuthal 
anisotropy of hadron spectra at high py, etc. [28]. On the other hand, one of 
the most spectacular features of low transverse momentum ( "soft" ) hadropro- 
duction at RHIC are strong collective flow effects: the radial flow (ordering of 
the mean transverse momentum of hadron species with corresponding mass) 
and the elliptic flow (mass-ordered azimuthal anisotropy of particle yields with 
respect to the reaction plane in non-central collisions). The development of 
such a strong flow is well described by the hydrodynamic models and requires 
short time scale and large pressure gradients, attributed to strongly interact- 
ing systems [29]. Note however, that results of hydrodynamic models usually 
disagree with the data on femtoscopic momentum correlations, resulting from 
the space-time characteristics of the system at freeze-out stage. 

The heavy ion collision energy in Large Hadron Collider (LHC) at CERN 
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a factor of 30 larger then that in RHIC, thereby allows one to probe new 
frontiers of super-high temperature and (almost) net-baryon free QCD [8- 
11]. The emphasis of the LHC heavy ion data analysis (at \/s = 5.5 TeV 
per nucleon pair for lead beams) will be on the perturbative, or hard probes 
of the QGP (quarkonia, jets, photons, high-p T hadrons) as well as on the 
global event properties, or soft probes (collective radial and elliptic flow effects, 
hadron multiplicity, transverse energy densities and femtoscopic momentum 
correlations). It is expected that at LHC energies the role of hard and semi- 
hard particle production will be significant even for the bulk properties of 
created matter. 

Another domain of QCD phase diagram is high baryon density region, which 
can be probed in the future experimental studies at relatively moderate heavy 
ion beam energies y/s ~ 10 GeV per nucleon pair (projects CBM [30] and 
MPD [31] at accelerator facilities FAIR-GSI and NICA-JINR, programs for 
heavy ion runs with lower beam energies at SPS and RHIC). These stud- 
ies are motivated by the intention to search for the "critical point" of the 
quark-hadron phase transition, predicted by lattice QCD, where the type of 
transition is changing from smooth "crossover" (at high temperatures and low 
net-baryon densities) to first order transition (at high net-baryon densities 
and low temperatures) [32]. 

Ongoing and future experimental studies of relativistic heavy ion collisions 
in a wide range of beam energies require the development of new Monte- 
Carlo (MC) event generators and improvement of existing ones. Especially for 
experiments which will be conducted at LHC, due to very high parton and 
hadron multiplicities, one needs fast (but realistic) MC tools for heavy ion 
event simulation [5-7] . The main advantage of MC technique for the simula- 
tion of high-multiplicity hadroproduction is that it allows a visual comparison 
of theory and data, including if necessary the detailed detector acceptances, 
responses and resolutions. A realistic MC event generator should include a 
maximum possible number of observable physical effects which are important 
to determine the event topology: from the bulk properties of soft hadropro- 
duction (domain of low transverse momenta pr ^ lGeV/c) such as collective 
flows, to hard multi-parton production in hot and dense QCD-matter, which 
reveals itself in the spectra of high-pr particles and hadronic jets. 

In most of the available MC heavy ion event generators, the simultaneous 
treatment of collective flow effects for soft hadroproduction and hard multi- 
parton in-medium production is absent. For example, the popular MC model 
HIJING [33] includes jet production and jet quenching on some level, but 
it does not include any significant flow effects. The recently developed MC 
model JEWEL [34], which combines a perturbative final state parton shower 
with QCD-medium effects, is actually not a heavy ion event generator, but 
rather a simulator of jet quenching in individual nucleon- nucleon sub-collisions 
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(like PYQUEN [12,14]). The event generators FRITIOF [35] and LUCIAE [36] 
include jet production (but without jet quenching), while some collective nu- 
clear effects (such as string interactions and hadron rescatterings) are taken 
into account in LUCIAE. Another MC model THERMINATOR [37] includes 
detailed statistical description of "thermal" soft particle production and can 
reproduce the main bulk features of hadron spectra at RHIC (in partic- 
ular, describe simultaneously the momentum-space measurements and the 
freeze-out coordinate-space data), but it does not include hard parton pro- 
duction processes. There is a number of microscopic transport hadron models 
(UrQMD [38], QGSM [39], AMPT [40], etc.), which attempt to analyze the soft 
particle production in a wide energy range, however they also do not include 
in-medium production of high-j>r multi-parton states. Moreover, hadronic cas- 
cade models have difficulties with the detailed description of relatively low-pr 
RHIC data on elliptic flow and the size of hadron emission, obtained from 
momentum correlation measurements (e.g. AMPT can reproduce the elliptic 
flow or the correlation radii using different sets of model parameters). Another 
heavy ion event generator, ZPC [41], has been created to simulate parton cas- 
cade evolution in ultrarelativistic heavy ion collisions. From a physical point 
of view, such an approach seems reasonable only for very high beam energies 
(RHIC, LHC). Note also that the full treatment of parton cascades requires 
enormous amount of CPU run time (especially for LHC). 

Thus, in order to analyze existing data on low and high-pr hadron production, 
test the sensitivity of physical observables at the upcoming LHC experiments 
(and other future heavy ion facilities) to the QCD-matter formation, and study 
the experimental capabilities of constructed detectors, the development of ad- 
equate and fast MC models for simultaneous collective flow and jet quenching 
simulation is necessary. HYDJET++ event generator includes detailed treat- 
ment of soft hadroproduction as well as hard multi-parton production, and 
takes into account medium- induced parton rescattering and energy loss. Al- 
though the model is optimized for very high energies of colliding nuclei (RHIC, 
LHC), in practice it can also be used for studying the multi-particle produc- 
tion at lower energies at other heavy ion experimental facilities as FAIR and 
NICA. Moving from very high to moderately high energies, the contribution 
of the hard part of the event becames smaller, while the soft part turns into 
just a multi-parameter fit to the data. The heavy ion event in HYDJET++ 
is the superposition of two independent parts: the soft, hydro- type state and 
the hard state resulting from multi-parton fragmentation. Note that a con- 
ceptually similar approximation has been developed in [42] but, to our best 
knowledge, it is not implemented as an MC event generator. 

The main program of HYDJET++ is written in the object-oriented C++ 
language under the ROOT environment [1]. The hard part of HYDJET++ is 
identical to the hard part of Fortran-written HYDJET (version 1.5) and it is 
included in the generator structure as a separate directory. The soft part of 
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HYDJET++ event is the "thermal" hadronic state generated on the chemi- 
cal and thermal freeze-out hypersurfaces obtained from a parameterization of 
relativistic hydrodynamics with preset freeze-out conditions. It includes the 
longitudinal, radial and elliptic flow effects and the decays of hadronic reso- 
nances. The corresponding fast MC simulation procedure based on C++ code, 
FAST MC [17,18], was adapted to HYDJET++. The high generation speed 
in FAST MC is achieved due to almost 100% generation efficiency because 
of nearly uniform residual invariant weights which appear in the freeze-out 
momentum and coordinate simulation. 

Let us indicate some physical restrictions of the model. HYDJET++ is only 
applicable for symmetric AA collisions of heavy (A > 40) ions at high energies 
(V$ ~ 10 GeV). Since the hydro-type approximation for heavy ion collisions 
is considered to be valid for central and semi-central collisions, the results 
obtained for very peripheral collisions (with impact parameter of the order 
of two nucleus radii, b ~ ^Ra) m &y be not adequate. Nor do we expect a 
correct event description in the region of very forward rapidities, where the 
other mechanisms of particle production, apart from hydro-flow and jets, may 
be important. 



2 Physics model 

A heavy ion event in HYDJET++ is a superposition of the soft, hydro-type 
state and the hard state resulting from multi-parton fragmentation. Both 
states are treated independently. 



2. 1 Model for the hard multi-jet production 

The model for the hard multi-parton part of HYDJET++ event is the same 
as that for H YD JET event generator. A detailed description of the physics 
framework of this model can be found in the corresponding paper [12]. The 
approach to the description of multiple scattering of hard partons in the dense 
QCD-matter (such as quark-gluon plasma) is based on the accumulative en- 
ergy loss via the gluon radiation being associated with each parton scattering 
in the expanding quark-gluon fluid and includes the interference effect (QCD 
analog [43-47] of the well known Landau-Pomeranchuk-Migdal (LPM) effect 
in QED [48,49] for the emission of gluons with a finite formation time) using 
the modified radiation spectrum dE/dl as a function of decreasing tempera- 
ture T. The basic kinetic integral equation for the energy loss AE as a function 
of initial energy E and path length L has the form 
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where / is the current transverse coordinate of a parton, dPj dl is the scattering 
probability density, dE /dl is the energy loss per unit length, A = l/(ap) is the 
in-medium mean free path, p oc T 3 is the medium density at the temperature 
T, a is the integral cross section for the parton interaction in the medium. 
Both collisional and radiative energy loss are taken into account in the model. 

The partonic collisional energy loss due to elastic scatterings is treated in 
high-momentum transfer limit [50-52]: 



dE co1 1 tm r da 

ll =WTa J < 2 > 
where the dominant contribution to the differential scattering cross section is 



da gTrojg) E 2 12vr 

dt t 2 E 2 -mf s (33 - 2N f ) In (t/A 2 QCD ) { ) 

for the scattering of a hard parton with energy E and mass m p off the 
"thermal" parton with energy (or effective mass) m ~ 3T <C E. Here 
C = 9/4, 1,4/9 for gg, gq and qq scatterings respectively, a s is the QCD run- 
ning coupling constant for Nf active quark flavors, and Aq C d is the QCD scale 
parameter which is of the order of the critical temperature of quark-hadron 
phase transition, Aqcd — T c ~ 200 MeV. The integrated cross section a is 
regularized by the Debye screening mass squared p 2 D {T) ~ 4ira s T 2 (l +Nf/6). 
The maximum momentum transfer t max = [s — {m p + mo) 2 ][s — {m p — m ) 2 ]/s 
where s = 2m -E + + m 2 . The model simplification is that the collisional 
energy loss due to scatterings with low momentum transfer (resulting mainly 
from the interactions with plasma collective modes or colour background fields, 
see for recent developments [53-61] and references therein) is not considered 
here. Note, however, that in the majority of estimations, the latter process 
does not contribute much to the total collisional loss in comparison with the 
high-momentum scattering (due to absence of the large factor ~ In (E /[/,&)), 
and in numerical computations it can be effectively "absorbed" by means of 
redefinition of minimum momentum transfer t m i n ~ p? D in (2). Another model 
assumption is that the collisional loss represents the incoherent sum over all 
scatterings, although it was argued recently in [62] that interference effects 
may appear in elastic parton energy loss in a medium of finite size. 

The partonic radiative energy loss is treated in the frameworks of BDMS 
formalism [63,64]. In fact, there are several calculations of the inclusive en- 
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ergy distribution of medium-induced gluon radiation using Feyman multiple 
scattering diagrams. The detailed discussions on the relation between these 
approaches, their basic parameters and physics predictions can be found in the 
recent proceedings of CERN TH Workshops (see [8,9] and references therein). 
In the BDMS approach, the strength of multiple scattering is characterized 
by the transport coefficient q = fi 2 D /\ g (X g is the gluon mean free path), 
which is related to the elastic scattering cross section a in (3). In our simu- 
lations this strength is regulated mainly by the initial QGP temperature To- 
Then the energy spectrum of coherent medium-induced gluon radiation and 
the corresponding dominant part of radiative energy loss of massless parton 
become [63,64]: 



E 




j du 


y 2 ] 

l-y + — 
y 2 



— = / du l-y + — In cos(wiTi) , (4) 

dl irL 



^f^-V+^V 2 )**™ with -=^f^y (5) 

where T\ = Lj (2A 9 ), y = uj/E is the fraction of the hard parton energy carried 
away by the radiated gluon, and Cr = 4/3 is the quark color factor. A similar 
expression for the gluon jet can be obtained by setting Cr = 3 and properly 
changing the factor in the square brackets in (4) [63]. The integration (4) is 
carried out over all energies from u m i n = -Elpm = H 2 d\, the minimum radiated 
gluon energy in the coherent LPM regime, up to initial parton energy E. Note 
that we do not consider here possible effects of double parton scattering [65,66] 
and thermal gluon absorption [67], which can be included in the model in the 
future. 

The simplest generalization of the formula (4) for a heavy quark of mass m q 
can be done by using the "dead-cone" approximation [68]: 



dE _ 1 dE 

dldZ^ ~ (1 + (Puff 2 ) 2 dldJ mq=Q ' 




(6) 



One should mention a number of more recent developments in heavy quark en- 
ergy loss calculations available in the literature (see [8] and references therein) , 
which can be included in the model in the future. 

The medium where partonic rescattering occurs is treated as a boost-invariant 
longitudinally expanding quark-gluon fluid, and the partons as being produced 
on a hyper-surface of equal proper times r [69]. In order to simplify numerical 
calculations we omit here the transverse expansion and viscosity of the fluid 
using the well-known scaling solution obtained by Bjorken [69] for a temper- 
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ature T and energy density e of QGP at T > T c ~ 200 MeV: 



e(r)r^ = e r 4/3 , T{t)t^ = T r 1/3 , p(r)r = p r . (7) 

The proper time tq of QGP formation, the initial QGP temperature To at mid- 
rapidity (y = 0) for central (impact parameter 6 = 0) Pb+Pb collisions, and 
the number Nf of active flavours are input model parameters. For non-central 
collisions and for other beam atomic numbers, T is calculated automatically: 
the initial energy density Eo{b) is supposed to be proportional to the ratio of 
nuclear overlap function T AA {b) to effective transverse area of nuclear overlap- 
ping S AA {b) [52]: 



£o(6) = 50(6 = 0) r AA (b = o) s AA (b) ' £o(A) = £o(Pb) Km) ' (8) 

T AA (b) = J J rdrT A ( ri )T A (r 2 ) , S AA (b) = j '# J rdr , (9) 





T A (r) = A J p A (r,z)dz , n j2 = y r 2 + j ± rbcostp , (10) 
i? cfr = minlyi?^ _ _ s i n 2 ^ + -cosip,]lR\ - — sin 2 ^ - -cos^}, (11) 

where r 1)2 (6, r, ^) are the distances between the centres of colliding nuclei 
and the jet production vertex V(r cosip, r sin^) (r being the distance from 
the nuclear collision axis to V), R c s{b,ip) is the transverse distance from the 
nuclear collision axis to the effective boundary of nuclear overlapping area in 
the given azimuthal direction ip (so that R c a(b = 0) is equal to the nuclear 
radius R A = 1.15A 1 / 3 ), T A (b) is the nuclear thickness functions, and p A (r,z) 
is the standard Woods-Saxon nucleon density distribution [70] (see [52] for 
detailed nuclear geometry explanations). The rapidity dependent spreading of 
the initial energy density around mid-rapidity y = is taken in the Gaussian- 
like form. The radial profile of the energy density (or the temperature) is also 
introduced in the parton rescattering model. The transverse energy density in 
each point inside the nuclear overlapping zone is proportional to the impact- 
parameter dependent product of two nuclear thickness functions T A in this 
point: 



e(n, r 2 ) cx T A (n)T A (r 2 ) . (12) 

Thus, in fact, the input parameter T oc e^ 4 acquires the meaning of an 
effective (i.e. averaged over the whole nuclear overlapping region) initial tem- 
perature in central Pb+Pb collisions. 

Note that the other scenarios of QGP space-time evolution for the MC imple- 
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mentation of the model were also studied. In particular, the influence of the 
transverse flow, as well as that of the mixed phase at T = T c , on the inten- 
sity of jet rescattering (which is a strongly increasing function of T) has been 
found to be inessential for high initial temperatures T 3> T c . On the contrary, 
the presence of QGP viscosity slows down the cooling rate, that implies a jet 
parton spending more time in the hottest regions of the medium. As a result 
the rescattering intensity increases, i.e., in fact the effective temperature of 
the medium appears to be higher as compared with the perfect QGP case. We 
also do not take into account here the probability of jet rescattering in the 
nuclear matter, because the intensity of this process and the corresponding 
contribution to the total energy loss are not significant due to much smaller 
energy density in a "cold" nucleus. 

Another important element of the model is the angular spectrum of in-medium 
gluon radiation. Since the detailed calculation of the angular spectrum of emit- 
ted gluons is rather sophisticated, the simple "small-angle" parameterization 
of the gluon distribution over the emission angle 9 was used: 



where 9q ~ 5° is the typical angle of the coherent gluon radiation as esti- 
mated in [71]. Two other parameterizations ("wide-angle" dN 9 /d9 oc 1/9, and 
"collinear" dN 9 /d9 = 5(9)) are also envisaged. 

The model for single hard nucleon-nucleon sub-collision PYQUEN [14] was 
constructed as a modification of the jet event obtained with the generator of 
hadron-hadron interactions PYTHIA_6.4 [15]. The event-by-event simulation 
procedure in PYQUEN includes generation of the initial parton spectra with 
PYTHIA and production vertexes at given impact parameter, rescattering-by- 
rescattering simulation of the parton path in a dense zone, radiative and col- 
lisional energy loss per rescattering, final hadronization with the Lund string 
model for hard partons and in-medium emitted gluons. Then the PYQUEN 
multi-jets generated according to the binomial distribution are included in the 
hard part of the event. If there is no nuclear shadowing, the mean number of 
jets produced in AA events at a given impact parameter b is proportional to 
the mean number of binary NN sub-collisions, N^n = TAAip)a l ^ N (y/s). In the 
presence of shadowing, this number is determined as 




(13) 




T A (ri)T A (r 2 )S(r u r 2 ,p T ,y) , 



(14) 
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where ctnn{V^) an d d&mf(pT, y/s) / ' dp\dy , calculated with PYTHIA, are the 
total inelastic non-diffractive NN cross section and the differential cross section 
of the corresponding hard process in NN collisions (at the same c.m.s. energy, 
y/s, of colliding beams) with the minimum transverse momentum transfer p™ m 
respectively. The latter is another input parameter of the model. In the HY- 
DJET frameworks, partons produced in (semi)hard processes with the mo- 
mentum transfer lower than p™ m are considered as being "thermalized" , so 
their hadronization products are included in the soft part of the event "auto- 
matically". The factor S < 1 in (14) takes into account the effect of nuclear 
shadowing on parton distribution functions. It can be written as a product of 
shadowing factors for both the colliding nuclei as 

S(n, r 2 ,p T , y) = S\( Xl , Q\ n) S{(x 2 , Q\ r 2 ) , (15) 

where S l A is the ratio of nuclear to nucleon parton distribution functions (nor- 
malized by the atomic number) for the parton of type (light quark or 
gluon), Xi t 2 are the momentum fractions of the initial partons from the incom- 
ing nuclei which participate in the hard scattering characterized by the scale 
Q 2 = £ix 2 s, and r 12 are the transverse coordinates of the partons in their 
respective nuclei, so that r x + r 2 = b. In HYDJET++, nuclear shadowing 
corrections are implemented not by modifying the parton showering of single 
NN collisions in PYTHIA but by correcting of the contribution of initial coher- 
ent, multiple scattering in effective way. In fact, this nuclear effect reduces the 
number of partons in the incoming hadronic wave-function of both the nuclei 
and thus reduces the total jet production cross section. Since the degree of 
this reduction depends on the kinematic variables of incoming hard partons, 
both initial and final parton momentum spectra are modified as a result of 
nuclear shadowing. 

The shadowing corrections are expected to be very significant at LHC ener- 
gies, where both soft and relatively high-p T (< 15 GeV/c) particle production 
probe the \ow-x gluon distribution of the target at moderate scales, Q 2 ~ p\, 
and are therefore strongly influenced by unitarity effects. In the Glauber- 
Gribov theory [72], this phenomenon arises from coherent interaction of the 
projectile fluctuation on the target constituents and is closely related to the 
diffractive structure function of the nucleon. Due to the factorization theorem 
for hard processes in QCD, S\ describes the modifications of nuclear parton 
distribution functions (i.e. distributions of quarks and gluons in nuclei), such 
that 

f i/A (x,Q 2 ,b) = f i/N (x,Q 2 )S t (A,b,x,Q 2 ) . (16) 

From summation of Pomeron fan diagrams the shadowing factor is found to 
be S\ = 1/ (1 + F l {x, Q 2 )T A (b)), where the effective cross section for quarks 
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and gluon, respectively, is found to be 



F\x, Q 2 ) = 4tt 



rO.l 

/ dxjp n 

J X 



Xjp < 



/^(/?,Q 2 )/£(:r,Q 2 ) 

(3g v ((3,Q 2 )/g(x,Q 2 ) 



(17) 
(18) 



where TP and g v denote the quark-singlet and gluon diffractive parton dis- 
tribution functions, £ and g are the normal parton distribution functions, B 
and fjp are the slope of diffractive distribution and the Pomeron flux factor 
respectively, Fa is the nuclear form factor. The quark and gluon diffractive 
distributions are taken from the most recent experimental parameterizations 
by the HI Collaboration [73], and the resulting shadowing factors calculated 
in [16] are implemented in HYDJET and HYDJET++. 



2.2 Model for the soft "thermal" hadron production 



The soft part of HYD JET++ event is the "thermal" hadronic state generated 
on the chemical and thermal freeze-out hypersurfaces obtained from a parame- 
terization of relativistic hydrodynamics with preset freeze-out conditions. The 
detailed description of the physics frameworks of this model can be found in 
the corresponding papers [17,18]. It is supposed that a hydrodynamic expan- 
sion of the fireball ends by a sudden system breakup at given temperature T ch 
and chemical potentials ubi^SiPq (baryon number, strangeness and electric 
charge respectively). In this case, the momentum distribution of the produced 
hadrons retains the thermal character of the (partially) equilibrated Lorentz 
invariant distribution function in the fluid element rest frame [74,75]: 

/fV°; T ch , 7- ) = — t 9i — — — , (19) 

ls ! exp([p*°-/i J ]/T ch )±l 

where p*° is the hadron energy in the fluid element rest frame, gi = 2Jj + 1 
is the spin degeneracy factor, 7 S < 1 is the (optional) strangeness suppression 
factor, nf is the number of strange quarks and antiquarks in a hadron i. The 
signs ± in the denominator account for the quantum statistics of a fermion 
or a boson, respectively. Then the particle number density Pi Q (T, fii) can be 
represented in the form of a fast converging series [17], that is, 

pTiT,^) = §- 2 m 2 T±^eM^)K 2 (^) , (20) 

where K 2 is the modified Bessel function of the second order and mi is the 
particle mass. 
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The mean multiplicity Ni of a hadron species i crossing the space-like freeze- 
out hypersurface cr(x) in Minkowski space is computed using effective thermal 
volume approximation [76,77]: 

= p?(T, / d 3 a^x)u^(x) = p- q (T, ^)V eS , (21) 

where the four-vector d 3 a fJi (x) = n^(x)d 3 a(x) is the element of the freeze-out 
hypersurface directed along the hypersurface normal unit four- vector n^{x) 
with a positively defined zero component (n°(x) > 0), d 3 a(x) = \d 3 a IJi d 3 a^\ 1 ^' 2 
is the invariant measure of this element. The probability that the produced 
grand canonical ensemble consists of particles is thus given by Poisson 
distribution around its mean value Nf 

- (N-) Ni 

P(N l ) = e W (-N l ) [ -^± r . (22) 

The chemical potential \X{ for any particle species % at the chemical freeze-out 
is entirely determined by chemical potentials JT q per unit charge, i.e., per unit 
baryon number B, strangeness S, electric charge (isospin) Q, charm C, etc. It 
can be expressed as a scalar product, /i« = qljl, where ql = {B iy Si, Qi, Ci, ...} 
and Ji = {JIb, fts, ]j>q, Jic, ■■}■ Assuming constant temperature and chemical po- 
tentials on the chemical freeze-out hypersurface, the total quantum numbers 
q = {B, S, Q, C, ...} of the thermal part of produced hadronic system with cor- 
responding V e s can be calculated as q = V e sJ2i p7*Qi- Thus, the potentials p, q 
are not independent, so taking into account baryon, strangeness and electrical 
charges only and fixing the total strangeness S and the total electric charge 
Q, Jig and JIq can be expressed through baryonic potential JIb- Therefore, the 
mean number of each particle and resonance species at chemical freeze-out is 
determined solely by the temperature T and the baryonic chemical potential 
JIb, which can be related within the thermal statistical approach using the 
following parameterization [78]: 

T(fi B ) = a-bJl 2 B - cjx% , Ji B {Vs NN ) = - ' — , (23) 

1 + ey/s NN 

where a = 0.166 ± 0.002 GeV, b = 0.139 ± 0.016 GeV" 1 , c = 0.053 ± 
0.021 GeV" 3 and d = 1.308 ± 0.028 GeV, e = 0.273 ± 0.008 GeV" 1 . 

Since the particle densities at the chemical freeze-out stage may be too high 
to consider the particles as free streaming (see, e.g., [79]), the assumption of 
the common chemical and thermal freeze-outs can hardly be justified, so a 
more complicated scenario with different chemical and thermal freeze-outs is 
implemented in HYDJET++ (T ch > T th ). Within the concept of chemically 
frozen evolution, particle numbers are assumed to be conserved except for 
corrections due to decay of some part of the short-lived resonances that can 
be estimated from the assumed chemical to thermal freeze-out evolution time. 
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Then to determine the particle densities p^ q (T th , pf 1 ) at the temperature of 
thermal freeze-out T th , the conservation of the particle fractions from the 
chemical to thermal freeze-out evolution time is supposed, and an additional 
input model parameter, effective pion chemical potential /i° ff th at thermal 
freeze-out, is introduced: 

pr(T c v,) p-wn 

P^\pf) &(T*,vif th ) • 1 } 

Assuming for the particles heavier than pions the Boltzmann approximation 
in (20), one deduces from (20) and (24) the chemical potentials of hadrons at 
thermal freeze-out: 

„* _ T th ln ( pT(T c \p?) p^\pf^) \ 

Note that it can no longer be expressed in the form pi = cfip, which is valid 
only for chemically equilibrating systems. 

At relativistic energies, because of the dominant longitudinal motion, it is con- 
venient to parameterize the fluid flow four-velocity {u°(x),u(x)} = j{l,v(x)} 
at a point x in terms of the longitudinal (z) and transverse (r±) fluid flow 
rapidities 

v u (x)- l \n l + Vz{x) p u (x)- 1 ln 1 + ^(^) cosh ^) m 
1 K J 2 l-v z {xY F y ' 2 1 - v ± {x) coshr] u {x) ' v ; 

where v± = \v±\ is the magnitude of the transverse component of the flow 
three-velocity v = {v±,v z } = {v± cos0 u , v± sin M , v z }, i.e., 

u^(x) = {cosh p u cosh i] u , sinh p u cos 0", sinh p u sin M , cosh p u sinh i] u } 
= {(1 + u 2 ± y/ 2 coshr] u ,u ± , (1 + u 2 ± ) 1/2 sinh r] u } , 

u± = 7-uj_ = cosh r] u ^±v± , 7j_ = coshp". The linear transverse rapidity profile 
is used and the momentum anisotropy parameter 5(b) is introduced here: 



u x 



= \Jl + 5(b) sinh p u cos , tt y = y'l — 5(b) sinh p u sin , 

*- = :s^' c " (i=0) ' < 28 > 

where is the spatial azimuthal angle of the fluid element, r is its radial 
coordinate, p" iax (b = 0) is the maximal transverse flow rapidity for central 
collisions, and Rf(b) is the mean-square radius of the hadron emission re- 
gion [18], which determines the fireball transverse radius R(b,(f)) in the given 
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azimuthal direction <p through the spatial anisotropy parameter e(b): 



Jl-e 2 (b) 

R{bA) = Rf(b) . v ' • (29) 

yjl + e(b) cos 20 

For this case, it is straightforward to derive the formula for the total effective 
volume of hadron emission from the hypersurface of proper time r=const [17,18] : 



2tt R{b,4>) r? max 

V cS = J d 3 a fl (t,x)u^t,x) = rjd ( f> f (n^)rdr J f(rj)drj : (30) 

cr(t,x) 



where (n M w M ) = cosh p u yl + 5(b) tanh 2 p u cos 2(f), and f(r]) is the longitudinal 
flow rapidity profile (Gaussian or uniform). 

The value V e Q is calculated in HYDJET++ only for central collisions (6 = 0), 
and then for non-central collisions it is supposed to be proportional to the 
mean number of nucleons-participants iVp art (6) (39): 



V^b)=V cS (b = 0) J^f ) . (31) 

iVpart(0 = 0) 

Since the fireball transverse radius Rf(b = 0) and the freeze-out proper time 
Tf(b = 0) (as well as its standard deviation A 17 (6 = 0) - the emission duration) 
are the input model parameters, it is straightforward to determine these values 
for non-central collisions through the effective volume V e s and input anisotropy 
parameters e and 5 at the given b: 




^<°>(iwf • a - ( ^ a - (o) (^)J (33) 

Note that such choice of centrality dependence for R/(b) and 17 (6) is inspired 
by the experimentally observed dependence of momentum correlation radii 
R S idc(b), R out (b), R long (b) oc iV^) 1 / 3 [80,81]. The detailed study of femto- 
scopic momentum correlations at RHIC and LHC (including their centrality 
dependence) is planned in the HYDJET++ frameworks for the future. 

The momentum and spatial anisotropy parameters 5(b) and e(b) can be treated 
independently for each centrality, or can be related to each other through the 
dependence of the elliptic flow coefficient f 2 (e, 5) (the second-order Fourier 
coefficient in the hadron distribution over the azimuthal angle f relatively to 
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the reaction plane) obtained in the hydrodynamical approach [82]: 



2(5 -e) , A , 

K (i-P)d-^ ■ < 34) 

Then using the predicted by hydrodynamical models (and observed at RHIC) 
proportionality of the value v 2 (b) and the initial ellipticity e (b) = bjlRj^ (i.e. 
v 2 oc eo), one can derrive the relation between 5(b) and e(b) (supposing also 
e oc e ): 

Jl + 4B(e + B) -1 
6=-^ , B = C(l-e 2 )e, e = ke , (35) 

where C and k are the independent on centrality coefficients, which should be 
specified instead of 5(b) and e(b) dependences. The best fit of RHIC data to 
the elliptic flow (Figure 6 in Section 6) is obtained with the values C = 2 and 
k = 0.175. 

The "thermal" hadronic state in HYDJET++ consists of stable hadrons and 
resonances produced from the SHARE particle data table [19], which contains 
360 particles (excluding any not well established resonance states). The decays 
of resonances are controlled by the decay lifetime 1/T, there T is the resonance 
width specified in the particle table, and these decays occur with the probabil- 
ity density Texp(— Tr) in the resonance rest frame. Then the decay products 
are boosted to the reference frame in which the freeze-out hypersurface was 
defined. The space-time coordinates of the decaying particle are shifted from 
its initial position on the decay length AtP/M (M and P are the decaying 
particle mass and four- momenta respectively). The branching ratios are also 
taken from the SHARE [19]. Only the two- and three-body decays are taken 
into account in the model. The cascade decays are also possible. Resonances 
are given the mass distribution according to a non-relativistic Breit-Wigner 

P(m)dm oc ] 2/A dm , ( 36 ) 

(m — m Q y + Am 2 / 4 

where mo and Am are the resonance nominal mass and width respectively. 
The Breit-Wigner shape is truncated symmetrically, | m — mo | < Am, with 
Am taken for each particle from the PYTHIA [15] (Am = for some narrow 
resonances which are not present in PYTHIA). 



3 Simulation procedure 



Before any event generation, HYDJET++ run starts from PYTHIA initializa- 
tion at the given c.m.s. energy per nucleon pair yjs (input parameter f SqrtS), 
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and then it calculates the total inelastic NN cross section cr^ N (y/s) (output 
parameter Sigin) and the hard scattering NN cross section cr^f^(\/s,p^ in ) 
(output parameter Sigjet) with the minimum transverse momentum trans- 
fer p™ m (input parameter fPtmin). Then the tabulation of nuclear thickness 
function T A (b) and nuclear overlap function T AA (b) is performed. If the impact 
parameter b of heavy ion AA collision is not fixed (input parameter f If b 7^ 0), 
its value b (output parameter Bgen) is generated in each event between the 
minimum (input parameter f Bmin) and maximum (input parameter fBmax) 
values in accordance with the differential inelastic AA cross section, obtained 
from the standard generalization of Glauber multiple scattering model [83] to 
the case of independent inelastic nucleon-nucleon collisions: 



d 2 (T': 



d 2 b 



<b, y/i) = 



1 - ^ 2 TAA(b)a^ N (V~s) ^ 



(37) 



If the impact parameter b is fixed (flfb = 0), then its value Bgen in each 
event is equal to the input parameter fBf ix. After specification of b for each 
given event, the mean numbers of binary NN sub-collisions N hin (output pa- 
rameter Nbcol) and nucleons-participants N part (output parameter Npart) are 
calculated: 



N U n(b, >/£) = T AA (b)a™ N (^) , 



2tt 



(38) 



AW(&, V~s)=Jd^jJ rdrT A ( ri ) [l - exp {o% N (VZ)T A (r 2 )}] . (39) 



The next step is the simulation of particle production proper in the event. 
The soft, hydro-type state and the hard, multi-jet state are simulated inde- 
pendently. When the generation of soft and hard states in each event at given 
b is completed, the event record (information about coordinates and momenta 
of initial particles, decay products of unstable and stable particles) is formed 
as the junction of these two independent event outputs in the RunOutput . root 
file. 



3.1 Generation of the hard multi-jet state 



The following event-by-event MC simulation procedure is applied to generate 
the hard state resulting from multi-parton fragmentation in the case when jet 
production is switched on (input parameter fNhsel=l, 2, 3 or 4). 

(1) Calculation of the number of NN sub-collisions N 3 ^ A (output parameter 
Njet) producing hard parton-parton scatterings of selected type (QCD- 
dijets by default, PYTHIA parameter msel=l) with p T > p™ m , according 
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to the binomial distribution around its mean value (14) (without shadow- 
ing correction yet, 5 = 1). For this purpose, each of Nbcol sub-collisions 
is tested basing on comparison of random number generated uniformly 
in the interval [0, 1] with the probability pjet=Sigjet/Sigin to produce 
the hard process. The i-th sub-collision is accepted if & <pjet, and is 
rejected in the opposite case. 

(2) Selecting the type of hard NN sub-collision (pp, np or nn) in accordance 
with the phenomenological formula for the number of protons Z in the 
stable nucleus A, Z = A/ (1.98+0. 015 A 2 / 3 ). For this purpose, each of Nj et 
"successful" sub-collisions is tested basing on comparison of two random 
numbers and generated uniformly in the interval [0, 1] with the prob- 
ability Z/A. The proton-proton sub-collision is selected if < Z/A, 
neutron-neutron sub-collision — if > Z/A, and proton- neutron sub- 
collision — in other cases. 

(3) Generation of multi-parton production in Njet hard NN sub-collisions 
by calling PYTHIA Njet times (call u pyinit and call u pyevnt, parton 
fragmentation being switched off by setting mstp(lll)=0). The spatial 
vertex of jet production for each sub-collision is generated by PYQUEN 
routine according to the distribution 

d^rdr {) T AA (b) { ' 



(4) If jet quenching is switched on (f Nhsel=2 or 4), initial PYTHIA-produced 
partonic state is modified by QCD-medium effects with PYQUEN gener- 
ator for each sub-collision separately. The radiative and collisional energy 
loss are both taken into account by default (input parameter f lenglu=0), 
but the options to have only radiative loss (f Ienglu=l) or only collisional 
loss (f Ienglu=2) are envisaged. For each hard parton with the initial 
transverse momentum p T > 3 GeV/c and pseudorapidity | rj |< 3.5 in 
the given sub-collision, the following rescattering scheme is applied. 

• Calculation of the scattering cross section crfa) — J dt da/dt and gen- 
eration of the transverse momentum transfer t(r«) in the i-th scattering 
according to (3) (Tj is a current proper time). 

• Generation of the displacement between the i-th and (i + l)-th scatter- 
ings, k = (r i+ i -Tt): 

dP h r 

— = exp (- / A- 1 ^ + s)ds) , A-V) = a(r)p(r) , (41) 

1 o 

and calculation of the corresponding transverse distance, I^pt/E. 

• Generation of the energy of a radiated in the i-th scattering gluon, 
uji = AE raxi: i, according to (4) and (6): 
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dJ, _2a a (ji 2 D )XC R \ 

V 2 



, |m„=U r 



In |cos(u;iTi)| , (42) 



dl 1 dl 

and its emission angle Oi relative to the parent parton determined accord- 
ing to the small-angle distribution (13) (input parameter flanglu=0). 
The options to have wide-angle (f Ianglu=l) or collinear (f Ianglu=2) 
distributions are also envisaged. 

• Calculation of the collisional energy loss in the i-th scattering 

AE„ W = , (44) 

where energy of "thermal" medium parton m is generated according to 
the isotropic Boltzmann distribution at the temperature T{ji). 

• Reducing the parton energy by collisional and radiative loss per each 
i-th. scattering, 

A£ toM = A£ coM + A£ radiJ , (45) 

and changing the parton momentum direction by means of adding the 
transverse momentum recoil due to elastic scattering i, 

Ak 2 tl = (E — -^-f - (p - - _ A)2 _ m 2 (46) 
1 2mJ ^ p2m 0t 2p> p 1 ; 

• Going to the next rescattering, or halting the rescattering if one of the 
following two conditions is fulfilled: (a) the parton escapes the hot QGP 
zone, i.e. the temperature in the next point T(r i+ i, r i+1 , t^+i) becomes 
lower than T c = 200 MeV; or (b ) the parton loses so much of energy that 
its transverse momentum pr^+i) drops below the average transverse 
momentum of the "thermal" constituents of the medium, 2T{ji + i). In 
the latter case, such a parton is considered to be "thermalized" and its 
momentum in the rest frame of the quark-gluon fluid is generated from 
the random "thermal" distribution, dN/d 3 p oc exp (— E/T), boosted to 
the center-of-mass of the nucleus-nucleus collision. 

At the end of each NN sub-collision, adding new (in-medium emitted) glu- 
ons to the PYTHIA parton list and rearrangement of partons to update 
string formation with the subroutine PYJOIN are performed. An addi- 
tional gluon is included in the same string as its "parent", and colour 
connections of such gluons are re-ordered relative to their ^-coordinates 
along the string. 

(5) If the nuclear shadowing is switched on (input parameter f Ishad=l) and 
beam ions are Pb, Au, Pd or Ca (fAw=207., 197., 110. or 40.), each 
hard NN sub-collision is tested basing on comparison of random number 
& generated uniformly in the interval [0, 1] with the shadowing factor S 
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(15) taken from the available parameterization (subroutine u ggshad). It 
is determined by the type of initially scattered hard partons, momen- 
tum fractions obtained by the partons at the initial hard interaction 
£1,2 (PYTHIA parameters pari (33), pari (34)), the square of trans- 
verse momentum transfer in the hard scattering Q 2 (PYTHIA parameter 
pari (22)), and the transverse position of jet production vertex relative 
to the centres of nuclei r\p- The given sub-collision is accepted if £j < S, 
and is rejected in the opposite case. 

(6) Formation of hadrons for each "accepted" hard NN sub-collision with 
PYTHIA (parton fragmentation being switched on by call u pyexec), 
and final junction of Njet sub-events to common array hyjets using 
standard PYTHIA event output format. If particle decays in PYTHIA 
are switched off (input PYTHIA parameter mst j (21)=0), then forma- 
tion of final hadrons from the two- and three-body decays of resonances 
follows the same pattern as that for soft hydro-part with the branching 
ratios taken from the SHARE particle decay table [19]. 



3.2 Generation of the soft "thermal" state 

The following event-by-event MC simulation procedure is applied to generate 
the soft "thermal" state in the case when soft hadroproduction is switched on 
(input parameter fNhsel=0, 1, or 2). 

(1) Initialization of the chemical freeze-out parameters. It includes the cal- 
culation of particle number densities according to (20) (the calculation 
of chemical freeze-out temperature, baryon potential and strangeness po- 
tential as a function of beam c.m.s. energy f SqrtS according to (23) is 
performed in advance if the corresponding flag fTMuType>0). So far, only 
the stable hadrons and resonances consisting of u, d, s quarks are taken 
into account from the SHARE particle data table [19]. 

(2) Initialization of the thermal freeze-out parameters (if T th < T ch ). It in- 
cludes the calculation of chemical potentials according to (25) and particle 
number densities according to (24). 

(3) Calculation of the effective volume of hadron emission region V e &{b) ac- 
cording to (31) and (30), the fireball transverse radius Rf(b), the freeze- 
out proper time Tf(b) and the emission duration At/(6) according to (32) 
and (33), and the mean multiplicity of each particle species according to 
(21). Then the multiplicity is generated around its mean value according 
to the Poisson distribution (24). 

(4) For each hadron the following procedure to generate its four-momentum 
is applied (resonance mass is generated according to (36)). 
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• Generation of four-coordinates of a hadron in the fireball rest frame 

= {r cosh 77, rcos0, rsin0, rsinh??} on each freeze-out hypersurface 
segment r(r) for the element d 3 a^u^ = d 3 a^ = n^r) | 1 — (dr/dr) 2 I 1 / 2 
r(r)d 2 rdr], assuming and r to be the functions of r (i.e., indepen- 
dent of 77, 0). It includes sampling uniformly distributed in the interval 
[0, 2tt], generating 77 according to the uniform distribution in the interval 
[— Imax, ?7max] (if hag f EtaType<0) or Gaussian distribution exp(—r] 2 /2r] 2 ^ ax ) 
(if flag fEtaType>0) and r in the interval [0,R f (b)}) using a 100% effi- 
cient procedure similar to the ROOT routine GetRandomO . 

• Calculation of the corresponding collective flow four-velocities accord- 
ing to (27) and (28). 

• Generation of the three-momentum of a hadron in the fluid element 
rest frame p*{ sin 9*cos(f>*, sin 9*sin(f)*, cos 9*} according to the equi- 
librium distribution function fP(p°*-,T, fi^p* 2 dp* d cos 6* d<p* by means of 
sampling uniformly distributed cos 9* in the interval [—1,1] and 0* in 
the interval [0, 2ir], and generating p* using a 100% efficient procedure 
(similar to ROOT routine GetRandomO). 

• The standard von Neumann rejection/ acceptance procedure to take into 
account the difference between the true probability W* i d 3 ad 3 p*/p°* and 
the probability n°* fP(p°*;T, ^i)d 2 fdrjd 3 p* corresponding to the previous 
simulation steps. For this purpose, the residual weight is calculated [17]: 

W*d 3 a 

W res = — = 

n°*p *fi q d 2 rdi] 

Then the simulated hadron four-coordinate and four-momentum is tested 
basing on comparison of W- es with the random number generated 
uniformly in the interval [0, max(lV i res )]. The i-th hadron is accepted if 
& < W- es , and rejected in the opposite case (then the generation of its 
four-coordinate and four-momentum is repeated). 

• Boost of the hadron four-momentum in the center mass frame of the 
event using the velocity field v(x), that is, 

p = 7 (p°* + ^*) , p = p* + 7 (l + 7 )- 1 (p°* +p °)v . (48) 

Note that the high generation speed for this algorithm is achieved due 
to almost 100% generation efficiency because of nearly uniform residual 
weights Wl cs (47). 

(5) Formation of final hadrons from the two- and three-body decays of reso- 
nances with the random choice of decay channel according to the branch- 
ing ratios taken from the particle data file tabledecay.txt (by default, 
if flags f Decay and fWeakDecay >0). The particle "decay" programs de- 
veloped by N.S. Amelin for FAST MC generator [17,18] are implemented 
in HYDJET++. 

• For the two-body decay, the momenta |_px,2 1 of decay products in the 



n*p* 



(47) 
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resonance rest frame is calculated as: 



| p lj2 |= 0.5y / ((M 2 -m\- ml) 2 - 4mfm|)/M , (49) 

where M is the resonance mass, and ra^ are the decay product masses. 
The space orientation of the \p\p\ is generated randomly by sampling 
uniformly distributed cosine of the polar angle cos (9) in the interval 
[—1, 1] and the azimuthal angle in the interval [0, 27r] (in the resonance 
rest frame). The opposite directions for pi i2 vectors are chosen according 
to momentum conservation law. 

• For the three-body decay, the resonance energy is divided between 
the kinetic energy of three decay products uniformly (assuming constant 
matrix element of the decay): 

Ef n = 6AM , | pi |= v /((E 1 kin ) 2 + 2Ef a m 1 ) 

Ef n = (1 - 6)AM , \p 2 |= v /((^ 2 kin ) 2 + 2 J B 2 kin m 2 ) 
^3 in =(6-6)AM, | p 3 h V / ((^3 kin ) 2 + 2i?3 kin ^3) , 

where £i and £ 2 are the two random numbers distributed uniformly in 
the interval [0, 1] under the condition £ 2 > £i, and AM = M — mi — m 2 — 
m 3 . For the first particle, the space orientation of the \pi\ is generated 
randomly by sampling uniformly distributed cos(#) in the interval [—1,1] 
and in the interval [0, 27r] (in the resonance rest frame). Then the three- 
momenta of remaining two particles are calculated taking into account 
complanarity of the decay and the conservation laws: 

P2 H Vi I (sin#N cos 4>n cos 9 cos — sin 9n sin 0at sin + 
cos#Arsin0cos0) , 

P2 —\ V2 I (sin 9m cos0jv cos 9 sin — sin 9n sin 0jv cos0 + 

cos9 N sin 9 sin 0) , 
P2 =| P2 I (— sin6'7vcos0Arsin6' + cos 9 N cos 0) , 

P3 = ~(P2 +Pl) , 

COS^at = (| P2 | 2 - I P3 | 2 - I Pi H/(2 I Pi || p 3 |) , 

and 0at is generated uniformly in the interval [0,27r]. Finally, the mo- 
menta of decay particles for two- and three-body decays are boosted to 
the center-of-mass of the nucleus-nucleus collision with the resonance ve- 
locity. The space-time coordinates of the decay products coincide with 
the coordinates of decay of the initial resonance. 
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particles.data, 
tabledecay.txt 

(particle properties) 




RunlnputHy dj et 

(input model parameters) 



RunOutput.root 

(particle output information for each 
event and global output parameters) 



Tree "td' 

(particle output) 



ROOT macros 

(to produce histograms) 



root; 

Histograms 



Fig. 1. The block structure of HydjetH — h 
Overview of HYDJET++ software structure 



The basic frameworks of HYDJET++ are preset by the object-oriented C++ 
language and the ROOT environment [1]. There is also the Fortran- written 
part [13] which is included in the generator structure as a separate directory. 
The block structure of HYDJET++ is shown in Figure 1. The main program 
elements (particle data files, input and output files, C++ and Fortran routines) 
are described in details in this section below. 



4-1 Particle data files 



The information regarding the particle species included in the soft part of the 
HYDJET++ event is stored in the files particles . data and tabledecay . txt. 
The particles.data file contains the definition (PDG code) and physical 
properties (mass, decay width, spin, isospin, valence quark composition) of 
360 stable hadrons and resonances. The tabledecay.txt file contains decay 
channels and branching ratios. The structure of these files is the same as that in 
SHARE particle data table [19] and in event generator THERMINATOR [37], 
where the corresponding description in more details can be found. 

The file particles.data has the following format. 
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name - the particle label; 
mass - mass in GeV/c 2 ; 
width - width in GeV/c 2 ; 
spin - spin; 
I - isospin; 

13 - third component of isospin; 

q - number of light valence quarks in the particle; 

s - number of strange valence quarks in the particle; 

aq - number of light valence antiquarks in the particle; 

as - number of strange valence antiquarks in the particle; 

c - number of charm valence quarks in the particle; 

ac - number of charm valence antiquarks in the particle; 

MC - the particle identification code. 

The file tabledecay.txt has the following format. 

PdgParent - parent particle code; 

PdgDaughterl - first particle (decay product) code; 

PdgDaughter2 - second particle (decay product) code; 

PdgDaughter3 - third particle (decay product) code (appears for the three- 
body decays only); 

BR - the branching ratio for the decay. 

Note that the default PYTHIA_6.4 particle data settings [15] are used to 
generate hard part of HYDJET+- 1- event. 

4-2 Input parameters and files 

Run of HYDJET++ is controlled by the file RunlnputHydj et for different 
type of input parameters. For current version of the generator, two additional 
files with the optimized parameters for Au+Au collisions at y/s = 200A GeV 
(RunlnputHydj etRHTC200, Fig. 2) and for Pb+Pb collisions at y/s = 550CL4 
GeV (RunlnputHydj etLHC5500, Fig. 3) are available. To use them as the 
input one, the user should change the name of the corresponding file to the 
RunlnputHydj et. The default parameters for RunlnputHydj etRHIC200 were 
obtained by fitting RHIC data to various physical observables (see Section 6 
for the details). The default parameters for RunlnputHydj etLHC5500 represent 
our rough extrapolation from RHIC to LHC energy, and of course they may 
be varied freely by the user. 

The following input parameters should be specified by the user. 

f Nevnt - number of events to generate; 

f SqrtS - beam c.m.s. energy per nucleon pair in GeV; 
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f Aw - atomic weight of nuclei; 

f Ifb - flag of type of centrality generation (=0: impact parameter is fixed 
by fBf ix value, ^0: impact parameter is generated in each event between 
minimum (f Bmin) and maximum (f Bmax) values according to the distribu- 
tion (37)); 

f Bmin - minimum impact parameter in units of nucleus radius Ra; 
f Bmax - maximum impact parameter in units of nucleus radius Ra; 
f Bf ix - fixed impact parameter in units of nucleus radius Ra- 

The following input parameters may be changed by user from their default 
values. 

f Seed - parameter to set the random number seed (=0: the current time is 
used to set the random generator seed, t^O: the value fSeed is used to set 
the random generator seed and then the state of random number generator in 
PYTHIA is MRPY(l)=fSeed) (default: 0). 

Parameters for soft hydro- type part of the event. 

f T - chemical freeze-out temperature T ch in GeV; 
f MuB - chemical baryon potential per unit charge \Tb in GeV; 
f MuS - chemical strangeness potential per unit charge \Ts i n GeV; 
f MuI3 - chemical isospin potential per unit charge \Tq in GeV; 
f TthFO - thermal freeze-out temperature T th in GeV; 

f Mu_th_pip - chemical potential of positively charged pions at thermal freeze- 
out fif th in GeV; 

fTau - proper time at thermal freeze-out for central collisions Tf(b = 0) in 
fm/c; 

f SigmaTau - duration of emission at thermal freeze-out for central collisions 
At/ (6 = 0) in fm/c; 

f R - maximal transverse radius at thermal freeze-out for central collisions 
R f (b = 0) in fm; 

f Ylmax - maximal longitudinal flow rapidity at thermal freeze-out 77 max ; 
fUmax - maximal transverse flow rapidity at thermal freeze-out for central 
collisions p™ ax (& = 0); 

f Delta - momentum azimuthal anisotropy parameter at thermal freeze-out 5 
(for given event centrality class); 

f Epsilon - coordinate azimuthal anisotropy parameter at thermal freeze-out 
e (for given event centrality class); 

f If DeltaEpsilon - flag to use calculated 5 and e (<0: specifed by user values 
for given event centrality class are taken, >0: calculated for given impact 
parameter according to the parameterization (35) values are used in each 
event) (default: 0); 

f Decay - flag to switch on/off hadron decays (<0: decays "off", >0: decays 
"on") (default: 0); 
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fWeakDecay - flag to switch on/off weak hadron decays (<0: decays "off", >0: 
decays "on") (default: 0); 

f EtaType - flag to specify longitudinal flow rapidity distribution (<0: uniform 
in the range [-f Ylmax , f Ylmax] , >0: Gaussian with the dispersion fYlmax) 
(default: 1); 

f TMuType - flag to use calculated chemical freeze-out temperature, baryon 
potential and strangeness potential as a function of f SqrtS (<0: specified by 
user values are taken, >0: calculated according to the parameterization (23) 
values are used) (default: 0); 

fCorrS - flag and value to include strangeness suppression factor j s with 
fCorrS value (0 <fCorrS< 1, if fCorrS< then it is calculated using its 
phenomenological dependence 7 S = 1 — 0.386 exp (— 1.23T ch /yU£) from [84]) 
(default: 1.). 

Parameters for treatment of hard multi-parton part of the event. 

f Nhsel - flag to switch on/off jet and hydro-state production (0: jet production 
"off" and hydro "on", 1: jet production "on" and jet quenching "off" and hydro 
"on", 2: jet production "on" and jet quenching "on" and hydro "on", 3: jet 
production "on" and jet quenching "off" and hydro "off" , 4: jet production 
"on" and jet quenching "on" and hydro" off") (default: 2); 
f I shad - flag to switch on/off impact parameter dependent nuclear shadowing 
for gluons and light sea quarks (u,d,s) (0: shadowing "off", 1: shadowing "on" 
for fAw=207., 197., 110. or 40.) (default: 1); 

f Ptmin - minimal transverse momentum transfer p™ in of hard parton-parton 
scatterings in GeV/c (the PYTHIA parameter ckin(3)=f Ptmin) 

PYQUEN energy loss model parameters: 

f TO - initial temperature T (in GeV) of QGP for central Pb+Pb collisions 

at mid-rapidity (initial temperature for other centralities and atomic numbers 

will be calculated automatically) (allowed range is 0.2<fT0<2.); 

fTauO - proper QGP formation time To in fm/c (0 . 01<f Tau0<10); 

fNf - number of active quark flavours Nf in QGP (fNf =0, 1, 2 or 3); 

f Ienglu - flag to fix type of in-medium partonic energy loss (0: radiative and 

collisional loss, 1: radiative loss only, 2: collisional loss only) (default: 0); 

f Ianglu - flag to fix type of angular distribution of in-medium emitted gluons 

(0: small-angular (13), 1: wide-angular, 2: collinear) (default: 0). 

Note that if specified by user value of input parameter extends out of the 
allowed range, its default value is used in HYDJET+- 1- run. 

A number of important PYTHIA parameters also may be changed/specified in 
RunlnputHydj et file (they are not shown in Figs. 2 and 3). The rest PYTHIA 
parameters can be changed (if it is necessary) using corresponding common 
blocks in the file progs_f ortran.f . 
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JVumber of events to generate 

C.m.s. energy per nucleon pair, fSqrtS [GeV] 
200. 

Atomic weigth of nuclei, fAw 
197. 

Flag of type of centrality generation, fBfix (=0 is fixed by fBfix, >0 distributed [iBfmin, fBmax]) 



Minimum impact parameter in units of nuclear radius, fBmin 



^Maximum impact parameter in units of nuclear radius, fBmax 
Fixed impact parameter in units of nuclear radius, fBfix 



Parameter to set random number seed, fSeed (=0 the current time is used, >0 the value fSeed is used) 



Temperature at chemical freeze-out, fT [GeV] 
0.165 

Chemical baryon potential per unit charge, fMuB [GeV] 

^Chemical strangeness potential per unit charge, fMuS [GeV] 

Chemical isospin potential per unit charge fMuI3, [GeV] 

Temperature at thermal freeze-out, fTthFO [GeV] 
0.1 

Chemical potential of pi+ at thermal freeze-out, fMu_th_pip [GeV] 
0.053 

Proper time proper at thermal freeze-out for central collisions, fTau [fm/c] 



Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c] 



Maximal transverse radius at thermal freeze-out for central collisions, fR [fm] 
10. 

Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax 
3.3 

Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax 
^Momentum azimuthal anizotropy parameter at thermal freeze-out, fDelta 

Spatial azimuthal anisotropy parameter at thermal freeze-out, fEpsilon 
^Flag to specify fDelta and fEpsilon values, flfDeltaEpsilon (<=0 user's ones, >0 calculated) 
^Flag to switch on/off hadron decays, fDecay (<0 decays off, >=0 decays on) 

Flag to switch on/off weak hadron decays, fWeakDecay: (<0 decays off, >=0 decays on) 

Flag to choose rapidity distribution, fEtaType (=<0 uniform, >0 Gaussian with the dispersion Ylmax) 



Flag to use calculated T_ch, mu_B and mu_S as a function of fSqrtS, fTMuType 
(<=0 user's ones, >0 calculated) 

Flag and value to include the strangeness supression factor gamma_s with fCorrS value 
(0<fCorrS <=1, if fCorrS <= 0. - it will be calculated) 

Flag to include jet (J)/jet quenching (JQ) and hydro (H) state production, fNhsel 
Flag to switch on/off nuclear shadowing, flshad (0 shadowing off, 1 shadowing on) 



^Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c] 

Initial QGP temperature for central Pb+Pb collisions in mid-rapidity, fTO [GeV] 

Proper QGP formation time in fm/c, fTauO (0.01<fTau0<10) 
0.4 

Number of active quark flavours in QGP, fNf (0, 1, 2 or 3) 

2 

Flag to fixe type of partonic energy loss, flenglu 
....... . jj. g . , , ...... 



(0 radiative and collisional loss, 1 radiative loss only, 2 collisional loss only) 
i 

Flag to fix type of angular distribution of in-medium emitted gluons, flanglu 
^0 small-angular, 1 wide-angular, 2 collinear). 

Fig. 2. The input parameter file RunInputHydjetRHIC200 (default). 
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^Jumber of events to generate 

C.m.s. energy per nucleon pair, fSqrtS [GeV] 

Atomic weigth of nuclei, f Aw 
207. 

Flag of type of centrality generation, fBfix (=0 is fixed by fBfix, >0 distributed [fBfmin, fBmax]) 



Minimum impact parameter in units of nuclear radius, fBmin 



^Maximum impact parameter in units of nuclear radius, fBmax 
Fixed impact parameter in units of nuclear radius, fBfix 



Parameter to set random number seed, fSeed (=0 the current time is used, >0 the value fSeed is used) 



Temperature at chemical freeze-out, fT [GeV] 
0.170 

Chemical baryon potential per unit charge, fMuB [GeV] 



^Chemical strangeness potential per unit charge, fMuS [GeV] 

^Chemical isospin potential per unit charge fMuI3, [GeV] 

Temperature at thermal freeze-out, fTthFO [GeV] 
0.13 

^Chemical potential of pi+ at thermal freeze-out, fMu_th_pip [GeV] 

Proper time proper at thermal freeze-out for central collisions, fTau [fm/c] 

Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c] 
3. 

^Maximal transverse radius at thermal freeze-out for central collisions, fR [fm] 

Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax 
4. 

^Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax 
^Momentum azimuthal anizotropy parameter at thermal freeze-out, fDelta 
^Spatial azimuthal anisotropy parameter at thermal freeze-out, fFpsilon 
Flag to specify fDelta and fEpsilon values, flfDeltaEpsilon (<=0 user's ones, >0 calculated) 
^Flag to switch on/off hadron decays, fDecay (<0 decays off, >=0 decays on) 
Flag to switch on/off weak hadron decays, fWeakDecay: (<0 decays off, >=0 decays on) 
^Flag to choose rapidity distribution, fEtaType (<= uniform, >0 Gaussian with the dispersion Ylmax) 





Flag to use calculated T_ch, mu_B and mu_S as a function of fSqrtS, fTMuType 
(<=0 user's ones, >0 calculated) 

Flag and value to include the strangeness supression factor gamma_s with fCorrS value 
(0<fCorrS <=1, if fCorrS <= 0., then it will be calculated) 

Flag to include jet (J)/jet quenching (JQ) and hydro (H) state production, fNhsel 
Flag to switch on/off nuclear shadowing, flshad (0 shadowing off, 1 shadowing on) 
Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c] 



Initial QGP temperature for central Pb+Pb collisions in mid-rapidity, fTO [GeV] 

^Proper QGP formation time in fm/c, fTauO (0.01<fTau0<10) 

^Number of active quark flavours in QGP, fNf (0, 1, 2 or 3) 

Flag to fixe type of partonic energy loss, flenglu 
illis" " " ' * " 



(0 radiative and collisional loss, 1 radiative loss only, 2 collisional loss only) 

Flag to fix type of angular distribution of in-medium emitted gluons, flanglu 
(0 small-angular, 1 wide-angular, 2 collinear). 



Fig. 3. The input parameter file RunInputHydjetLHC5500 (default). 
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4-3 Output parameters and files 



The program output is directed to a ROOT file specified by the user in the 
command line or by default to RunOutput .root. The output file contains a 
tree named td, which keeps the entire event record including primary particles 
and decay products with their coordinates and momenta information. Each 
decay product contains the unique index of its parent particle so that the 
entire event history may be obtained. Beside particle information, the output 
file contains also the following global output parameters for each event. 

Bgen - generated value of impact parameter b in units of nucleus radius Ra] 
Sigin - total inelastic NN cross section at given energy fSqrtS (in mb); 
Sigjet - hard scattering NN cross section at given fSqrtS, fPtmin (in mb); 
Ntot - generated value of total event multiplicity (Ntot=Nhyd+Npyt); 
Nhyd - generated multiplicity of "soft" hydro-induced particles; 
Npyt - generated multiplicity of "hard" jet-induced particles; 
N j et - generated number N 3 A e A of hard parton scatterings with p? >fPtmin; 
Nbcol - mean number of binary NN sub-collisions iVbi n (38) at given Bgen; 
Npart - mean number of nucleons-participants iVp art (39) at given Bgen. 

The event output tree (ROOT : : TTree) is organized as follows: 

td->Branch("nev" ,&nev, "nev/I") ; 

// event number 
td->Branch("Bgen" ,&Bgen, "Bgen/F") ; 

// generated impact parameter 
td->Branch( "Sigin" ,&Sigin, "Sigin/F") ; 

// total inelastic NN cross section 
td->Branch(" Sigjet" ,&Sigjet, "Sigjet/F") ; 

// hard scattering NN cross section 
td->Branch("Ntot" ,&Ntot, "Ntot/I") ; 

// total event multiplicity 
td->Branch("Nhyd" ,&Nhyd, "Nhyd/I") ; 

// multiplicity of hydro-induced particles 
td->Branch("Npyt" ,&Npyt , "Npyt/I") ; 

// multiplicity of jet-induced particles 
td->Branch("Njet" ,&Njet, "Njet/I") ; 

// number of hard parton-parton scatterings 
td->Branch( "Nbcol" , fcNbcol , "Nbcol/I " ) ; 

// mean number of NN sub-collisions 
td->Branch ( "Npart " ,&Npart, "Npart /I") ; 

// mean number of nucleon-participants 
td->Branch( "Px" , &Px [0] , "Px [npart] /F" ) ; 

// x-component of the momentum, in GeV/c 
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td->Branch ( "Py 11 , &Py [0] , "Py [npart] /F" ) ; 

// y-component of the momentum, in GeV/c 
td->Branch( "Pz" , &Pz [0] , "Pz [npart] /F" ) ; 

// z-component of the momentum, in GeV/c 
td->Branch("E",&E[0] , "E [npart] /F") ; 

// energy, in GeV 
td->Branch("X",&X[0] , "X [npart] /F") ; 

// x-coordinate at emission point, in fm 
td->Branch("Y",&Y[0] , "Y [npart] /F") ; 

// y-coordinate at emission point, in fm 
td->Branch("Z",&Z[0] , "Z [npart] /F") ; 

// z-coordinate at emission point, in fm 
td->Branch("T",&T[0] , "T [npart] /F") ; 

// proper time of particle emission, in fm/c 
td->Branch("pdg" ,&pdg[0] , "pdg [npart] /I") ; 

// Geant particle code 
td->Branch("Mpdg" ,&Mpdg[0] , "Mpdg [npart] /I") ; 

// Geant mother code (-1 for primordial) 
td->Branch("type" ,&type[0] , "type [npart] /I") ; 

// particle origin (=0: hydro, >0: jet) 
td->Branch ( " Index" , felndex [0] , " Index [Ntot] /I") ; 

// unique zero based index of the particle 
td->Branch( "Motherlndex" , Motherlndex [0] , "Motherlndex [Ntot] /I") ; 

// index of the mother particle (-1 if its a primary particle) 
td->Branch("NDaughters" ,&NDaughters [0] , "NDaughters [Ntot] /I") ; 

// number of daughter particles 
td->Branch( "Daughter llndex" , feDaughterllndex [0] , "Daughterl Index [Ntot] /I " ) 

// index of the first daughter (-1 if it does not exist) 
td->Branch( "Daughter2Index" , &Daughter2Index [0] , "Daughter2Index [Ntot] /I " ) 

// index of the second daughter 
td->Branch( "Daughter3Index" , &Daughter3Index [0] , "Daughter3Index [Ntot] /I " ) 

// index of the third daughter 

The possibility to create event output written directly to histograms (in ac- 
cording to user's specification in the file RunHadronSourceHisto . cxx) is envis- 
aged. The histogram output is directed by default to the file RunOutputHisto . root 
or to a file specified explicitly by the user in the command line. 



4-4 Organization of the Fortran package 

The hard, multi-jet part of HYDJET++ event is identical to the hard part 
of Fortran-written HYDJET [12,13] (version 1.5). The three Fortran files are 
placed in the directory PYQUEN. 
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pythia-6 .4. 11 . f - PYTHIA event generator to simulate hard NN sub-col- 
lisions [15] (is used if fNhsel=l, 2, 3 or 4). 

pyquenl_5 . f - PYQUEN event generator to modify PYTHIA-produced hard 
events according to the corresponding partonic energy loss model [12,14] (is 
used if f Nhsel=2 or 4), and to generate the spatial vertex of a jet production. 

progs_f ortran.f - main Fortran routine, which contains the following main 
subroutines and functions (the auxiliary subroutines and functions are not 
listed here). 

• subroutine u hyinit(fSqrtS,f Aw,f Ifb,f Bmin,fBmax,fBf ix) - initializes 
PYTHIA, calculates cross sections Sigin and Sigjet, tabulates T A (b), T AA (b). 

• subroutine u hyevnt - in each event generates the impact parameter value 
Bgen (if f Ifb ^ 0), calculates numbers Nbcol and Npart, generates number 
of "successful" hard NN sub-collisions Njet, and calls hyhard subroutine. 

• subroutine u hyhard - generates multi-parton production in Njet hard NN 
sub-collisions (with jet quenching if f Nhsel=2 or 4, and with nuclear shadow- 
ing if f Ishad=l), performs hadronization, and writes the final event output 
in the array hyjets. 

• subroutine u ggshad(inucl,x,q2,b,res,taf ) - adapts the parameteriza- 
tion of ratio of nuclear to nucleon parton distribution function res (res(l) 
for gluons and res (2) for sea quarks) in a given nucleus inucl, Bjorken x in 
the parton distribution, square of transverse momentum transfer in the hard 
scattering q2, and transverse position of jet production vertex relative to the 
nucleus center b (taf is the parametrized nuclear thickness function). 

• f unction u shadl (kf h,xbj ,q2,r) - in each event calculates the ratio of 
nuclear to nucleon parton distribution function (normalized by the atomic 
number) for the given parton code kfh, Bjorken x in a parton distribution 
xbj, square of transverse momentum transfer in the hard scattering q2, and 
transverse position of jet production vertex relative to the the nucleus center 
r (call u ggshad is performed). 

4-5 Organization of the C+ + package 

The organization of the C++ package in HYDJET++ is similar to the FAST 
MC framework presented in [17,18]. The C++ source files are placed in the 
main directory. 

The main modules of the C++ package. 
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• RunHadronSource . cxx, RunHadronSourceHisto . cxx - contains the main 
program. It handles the InitialStateHydjet class instance by calling the 
ReadParamsO, MultilniO, InitializeO andEvolveO member functions. 
The main program also creates the output ROOT tree and fills it with event 
information. 

• InitialState . cxx, InitialState . h - contains virtual class InitialState. 
The InitialState class is responsible for initialization of the PDG particle 
database (DatabasePDG). The function Evolve () organizes the particle decays 
and updates the coordinate-time information on the particle list. 

• InitialStateHydjet . cxx, InitialStateHydjet .h - contains the defini- 
tion of the class InitialStateHydjet. This class inherits the upper class 
InitialState. The member function ReadParamsO reads the parameters 
from the input file once per run. The function MultilniO calculates the total 
particle specie multiplicities once per run. The InitializeO function gener- 
ates the particles in the initial fireball. The Evolve () function, the heritage 
of the class InitialState, performs the decay of unstable particles. 

• HadronDecayer . cxx, HadronDecayer .h - contains two functions called by 
the function InitialState: : Evolve () to decay the resonances: Decay () and 
GetDecayTimeO. 

The service modules of the C++ package. 

• Particle . cxx, Particle. h - contains the Particle class which handles 
all the track information (PDG properties, momenta, coordinate, indexes of 
parent particles and decay products). 

• ParticlePDG . cxx, ParticlePDG.h - contains the definition of the class 
ParticlePDG which stores all the PDG properties of a particle specie. This 
class is used by the Particle class and by the DatabasePDG to organize the 
information. 

• DatabasePDG . cxx, DatabasePDG. h - contains the definition of the class 
DatabasePDG. This class reads and handles all the particle specie definitions 
found in particles . data and all the decay channels found in tabledecay . txt. 

• DecayChannel . cxx, DecayChannel .h - contains the DecayChannel class 
definition. This class stores the information for a single decay channel. A 
ParticlePDG instance uses DecayChannel objects to store informations re- 
garding decay channels. 

• GrandCanonical . cxx, GrandCanonical .h - contains the method to calcu- 
late the densities of energy, baryon and electric charge and particle numbers, 
within the grand canonical approach by means of the temperature and chem- 
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ical potentials. 

• StrangeDensity . cxx, StrangeDensity .h - contains the method to calcu- 
late the strangeness density within the grand canonical description. 

• HankelFunction. cxx, HankelFunction.h - computes modified Hankel func- 
tion of zero, first and second orders. 

• StrangePotential . cxx,StrangePotential . h - contains the method to cal- 
culate strange potential (if f TMuType u =l) from the initial strange density at 
given temperature and baryon potential. 

• EquationSolver . cxx, EquationSolver .h - is used by StrangePotential 

class to calculate strangeness potential. 

• RandArrayFunction. cxx, RandArrayFunction.h - defines several methods 
for shooting generally distributed random values, given a user-defined proba- 
bility distribution function. 

• UKUtility . cxx, UKUtility .h - contains the method IsotropicR3 to gen- 
erate the three-vector uniformily distributed on spherical surface, and the 
method MomAntiMom to calculate kinematical variables for two-body decays. 

• MathUtil.h - contains some useful constant determinations. 

• HY JET_C0MM0NS . h - contains the list of Fortran common blocks used by 
C++ package. 



5 Installation instructions 

HYDJET++ package is available via Internet. It can be downloaded from the 
corresponding web-page [85] as the archive HYDJET++ . ZIP, which contains the 
following files and directories: 

• the Makefile; 

• the input files (RunlnputHydjet, RunInputHydjetRHIC200, 
RunlnputHyd j etLHC5500) ; 

• HYDJET directory with inc (*.h files), src (*.cxx files) and SHARE data 
files (particle .data, tabledecay.txt); 

• PYQUEN directory with * . f files; 

• DOC directory (with short description and update notes for new versions); 

• RootMacros directory (the examples allowing one to obtain some physical 
results with HYDJET++, detailed comments are in the macros). 

In order to run HYDJET++ on Linux one needs: 
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1) C++ and Fortran compilers; 

2) ROOT libraries and include files. 

The main program is in the files RunHadronSource . cxx (for ROOT tree out- 
put) or RunHadronSourceHisto . cxx (for ROOT histogram output). To com- 
pile the package one needs to use the following commands in the main 
HYDJET directory: 
make (for ROOT tree output), or 

make u -f u Makef ile.HISTO (for ROOT histogram output). 
Then the executable file HYDJET (or HYDJET.HISTO) is created in the same 
directory. Once the program is compiled, one can use the executables above 
to start a simulation run. The input file necessary to run HYDJET++ must 
always be named RunlnputHydjet, so after preparing this file according to the 
examples included in the package one can start the simulation with the com- 
mand line: ./HYDJET (or . /HYDJET.HISTO). If an output file is not specified 
then the event records will be automatically directed to the file RunOutput . root 
(or RunOutputHisto . root). 



6 Validation of HYDJET++ with experimental RHIC data 



It was demonstrated in previous papers that FAST MC model can describe 
well the bulk properties of hadronic state created in Au+Au collisions at RHIC 
at y/s = 200A GeV (such as particle number ratios, low-p T spectra, elliptic 
flow coefficients v 2 (pt, b), femtoscopic correlations in central collisions) [17,18], 
while HYDJET model is capable of reproducing the main features of jet 
quenching pattern at RHIC (high-p^ hadron spectra and the suppression of az- 
imuthal back-to-back correlations) [12]. Since soft and hard hadronic states in 
HYDJET++ are simulated independently, a good description of hadroproduc- 
tion at RHIC in a wide kinematic range can be achieved, moreover a number 
of improvements in FAST MC and HYDJET have been done as compared to 
earlier versions. A number of input parameters of the model can be fixed from 
fitting the RHIC data to various physical observables (listed in Fig. 2). 

(1) Ratio of hadron abundances. It is well known that the particle abun- 
dances in heavy ion collisions in a wide energy range can be reasonable 
well described within statistical models based on the assumption that 
the produced hadronic matter reaches thermal and chemical equilibrium. 
The thermodynamical parameters fT B = 0.0285 GeV, fT s = 0.007 GeV, 
jiQ = —0.001, the strangeness suppression factor j s = 1, and the chemical 
freeze-out temperature T ch = 0.165 GeV have been fixed in [17] from fit- 
ting the RHIC data to various particle ratios near mid-rapidity in central 
Au+Au collisions at y/s = 200A GeV (7i-/tt + , p/n', K~/K + , K~ j-K', 
p/p, A/A, A/A, S/S, <j>/K~, A/p, S-/7T-). 
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(2) Low-pr hadron spectra. Transverse momentum pt and transverse mass 
itlt hadron spectra K + and p with m T < 0.7 GeV/c 2 ) near mid- 
rapidity at different centralities of Au+Au collisions at \/s = 200A GeV 
were analyzed in [18]. The slopes of these spectra allow the thermal freeze- 
out temperature T th = 0.1 GeV and the maximal transverse flow rapidity 
in central collisions p™ ax (& = 0) = 1.1 to be fixed. 

(3) Femtoscopic correlations. Because of the effects of quantum statistics 
and final state interactions, the momentum (HBT) correlation functions 
of two or more particles at small relative momenta in their c.m.s. are 
sensitive to the space-time characteristics of the production process on 
the level of fm. The space-time parameters of thermal freeze-out region 
in central Au+Au collisions at y/s = 200A GeV have been fixed in [18] 
by means of fitting the three-dimensional correlation functions measured 
for 7r + 7r + pairs and extracting the correlation radii -R S ide, -Rout and -Ri ong : 
Tf (b = 0) = 8 fm/c, Ar f {b = 0) = 2 fm/c, R f (b = 0) = 10 fm. 

(4) Pseudorapidity hadron spectra. The PHOBOS data on 77-spectra 
of charged hadrons [86] at different centralities of Au+Au collisions at 
yt~s = 200A GeV have been analyzed to fix the particle densities in 
the mid-rapidity region and the maximum longitudinal flow rapidity 
Vmax = 3.3 (Fig. 4). Since mean "soft" and "hard" hadron multiplicities 
depend on the centrality in different ways (they are roughly proportional 
to N part (b) and Abi n (6) respectively), the relative contribution of soft and 
hard parts to the total event multiplicity can be fixed through the central- 
ity dependence of dN/drj. The corresponding contributions from hydro- 
and jet-parts are determined by the input parameters fx'f th = 0.053 GeV 
and p™ m = 3.4 GeV/c respectively. 

(5) High-p T hadron spectra. High transverse momentum hadron spectra 
(jPt ^2 — 4 GeV/c) are sensitive to parton production and jet quenching 
effects. Thus fitting the measured high-p^ tail allows the extraction of 
PYQUEN energy loss model parameters. We assume the QGP formation 
time ro = 0.4 fm/c and the number of active quark flavours Nf = 2. 
Then the reasonable fit of STAR data on high-p^ spectra of charged pions 
at different centralities of Au+Au collisions at ^/s = 200A GeV [87] is 
obtained with the initial QGP temperature T = 0.3 GeV (Fig. 5). 

(6) Elliptic flow. The elliptic flow coefficient v 2 (which is determined as 
the second-order Fourier coefficient in the hadron distribution over the 
azimuthal angle <p relative to the reaction plane angle ipn, so that v 2 = 
(cos 2(tp — i/jr))) is an important signature of the physics dynamics at 
early stages of non-central heavy ion collisions. According to the typi- 
cal hydrodynamic scenario, the values v 2 (pt) at low-pr (~ 2 GeV/c) are 
determined mainly by the internal pressure gradients of an expanding 
fireball during the initial high density phase of the reaction (and it is 
sensitive to the momentum and azimuthal anisotropy parameters S and e 
in the frameworks of HYDJET++), while elliptic flow at high-p^ is gen- 
erated in HYDJET++ (as well as in other jet quenching models) due to 
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Fig. 4. The pseudorapidity distribution of charged hadrons in Au+Au collisions 
at y/s = 200A GeV for six centrality sets. The points are PHOBOS data [86], 
histograms are the HYD JET++ calculations (solid - total, dotted - jet part, dashed 
- hydro part). 
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Fig. 5. The transverse momentum distribution of positively charged pions in Au+Au 
collisions at y/s = 200A GeV for three centrality sets. The points are STAR 
data [87], histograms are the HYDJET++ calculations (solid - total, dotted - jet 
part, dashed - hydro part). 

the partonic energy loss in an azimuthally asymmetric volume of QGP. 
Figure 6 shows the measured by the STAR Collaboration transverse mo- 
mentum dependence of the elliptic flow coefficient v 2 of charged hadrons 
in Au+Au collisions at y/s = 200 A GeV for two centrality sets [88]. The 
values of 5 and e estimated for each centrality are written on the plots. 
Note that the choice of these parameters does not affect any azimuthally 
integrated physics observables (such as hadron multiplicities, r\- and pr- 
spectra, etc.), but only their differential azimuthal dependences. 
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Fig. 6. The transverse momentum dependence of the elliptic flow coefficient V2 of 
charged hadrons in Au+Au collisions at yfs = 200 A GeV for two centrality sets. 
The points are STAR data [88], histograms are the HYDJET++ calculations. 

7 Test run description 



The set of macros for comparison of various physics observables with the 
experimental data obtained for different centralities of Au+Au collisions at 
y/s = 200 A GeV are in the directory RootMacros. 

• f ig_eta_Phobos . C, f ig_eta_Phobos_read.C - the charged hadron pseudo- 
rapidity spectra are generated (Fig. 4). To start the generation, the user should 
type: root u -l u f ig_eta_Phobos . C, or root u -l u f ig_eta_Phobos_read.C. 

In the latter case, the number of events should be specified in the macro 
f ig_eta_Phobos_read.C. The PHOBOS data [86] are in the directory 
eta.PHOBOS. 

• f ig_PTH_STAR.C, f ig_PTH_STAR_read. C - the positive pion p^-spectra are 
generated (Fig. 5). To start the generation, the user should type: 
rootu-luf ig_PTH_STAR.C, or root u -l u f ig_PTH_STAR_read.C. In the latter 
case, the number of events should be specified in the macro 

f ig_PTH_STAR_read.C. The STAR data [87] are in the directory HPT.STAR. 

• f ig_v2_STAR. C, f ig_v2_STAR_read.C - the elliptic flow coefficients t>2(pr) 
of charged hadrons are generated (Fig. 6). To start the generation, the user 
should type: root u _ luf ig_v2_STAR. C, or root u _ luf ig_v2_STAR_read.C. In 
the latter case, the number of events should be specified in the macro 

f ig_v2_test_read.C. The STAR data [88] are in the directory v2_STAR. 
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8 Conclusion 



Ongoing and future experimental studies of relativistic heavy ion collisions re- 
quire the development of new Monte-Carlo event generators and improvement 
of existing ones. The main advantage of MC technique for the simulation of 
multiple hadroproduction is that it allows visual comparison of theory and 
data, including if necessary the detailed detector acceptances, responses and 
resolutions. The realistic MC event generator should include a maximum pos- 
sible number of observable physical effects which are important to determine 
the event topology: from bulk properties of soft hadroproduction (domain of 
low transverse momenta pr) to hard multi-parton production in hot QCD- 
matter, which reveals itself in spectra of high-pr particles and hadronic jets. 

The HYDJET++ event generator has been developed to simulate relativistic 
heavy ion A A collisions considered as a superposition of the soft, hydro- type 
state and the hard state resulting from multi-parton fragmentation. It in- 
cludes detailed treatment of soft hadroproduction (collective flow phenomena 
and resonance decays) as well as hard parton production, and takes into ac- 
count the known medium effects (as jet quenching and nuclear shadowing). 
The main program is written in the object-oriented C++ language under the 
ROOT environment. The hard part of HYDJET++ is identical to the hard 
part of Fortran-written HYDJET generator and it is included in the genera- 
tor structure as a separate directory. The soft part of HYDJET++ event is 
the "thermal" hadronic state generated on the chemical and thermal freeze- 
out hypersurfaces obtained from the parameterization of relativistic hydrody- 
namics with preset freeze-out conditions (the adapted C+ code FAST MC). 
HYDJET++ is capable of simultaneous reproducing the various hadronic ob- 
servables measured in heavy ion collisions at RHIC for different centrality sets 
in a wide kinematic range: ratio of hadron yields, pseudorapidity and trans- 
verse momentum spectra (for both low- and high-p^ domains), elliptic flow 
coefficients v 2 (pt), femtoscopic correlations. Although the HYDJET++ gen- 
erator is optimized for RHIC and LHC energies, in practice it can be also 
used for studying of multi-particle production in a wider energy range down 
to \/s ~ 10 GeV per nucleon pair at other heavy ion experimental facilities. 
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